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ABSTRACT 

Using only physical mechanisms, i.e., 3D MHD with no phenomenological 
viscosity, we have simulated the dynamics of a moderately thin accretion disk 
subject to torques whose radial scaling mimics those produced by lowest-order 
post-Newtonian gravitomagnetism. In this simulation, we have shown how, in the 
presence of MHD turbulence, a time-steady transition can be achieved between 
an inner disk region aligned with the equatorial plane of the central mass’s spin 
and an outer region orbiting in a different plane. The position of the equilibrium 
orientation transition is determined by a balance between gravitomagnetic torque 
and warp-induced inward mixing of misaligned angular momentum from the outer 
disk. If the mixing is interpreted in terms of diffusive transport, the implied 
diffusion coefficient is ~ (0.6-0.8)c 2 /D for sound speed c s and orbital frequency 
D. This calibration permits estimation of the orientation transition’s equilibrium 
location given the central mass, its spin parameter, and the disk’s surface density 
and scalehcight profiles. However, the alignment front overshoots before settling 
into an equilibrium, signaling that a diffusive model does not fully represent the 
time-dependent properties of alignment fronts under these conditions. Because 
the precessional torque on the disk at the alignment front is always comparable 
to the rate at which misaligned angular momentum is brought inward to the 
front by warp-driven radial motions, no break forms between the inner and outer 
portions of the disk in our simulation. Our results also raise questions about 
the applicability to MHD warped disks of the traditional distinction between 
“bending wave” and “diffusive” regimes. 
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Introduction 


It is now forty years since IBardeen fe Petterson (119751 1 pointed out that the precessional 
torques driven by the Lense-Thirring effect would, when applied to a fluid disk surrounding 
a spinning black hole, lead to alignment of its inner parts’ orbital axis with the black hole’s 
spin, while leaving untouched its original orientation, perhaps sharply misaligned, at large 
radius. Despite this long passage of time, there remains controversy over such basic issues as 
the detailed mechanism of alignment, where the transition between orientations occurs, and 
how its location depends on such disk_properties as the disk’s sound speed or surface density 
distribution. Bardeen fe Pettersonl (119751 1 suggested that the orientation transition would 
be determined by a competition between the Lense-Thirrin g torques and the imp ortati on of 
fresh misaligned angular momentum by the accretion flow. Papaloizou & Pringle (119831) rec¬ 
ognized that misaligned angular momentum can be more rapidly transported by the radial 
flows induced by disk-warping than by accretion, but the nature and speed of these radial 
flows was left uncertain because their regulation was ascribed to an imagined isotropic viscos¬ 
ity, which wa s assu me d to b e related to the local pressure through a universal dimensionless 
constant, a. [Pringle (119921) developed a simpler formalism to describe radial flow of mis¬ 
aligned angular momentum, in which it was supposed to spread diffusively, but the “diffusion 
coefficient” relating a radial gradient of angular momentum orientation to an associated an¬ 
gular momentum fl ux was e stimated assuming that radial flows are controlled by that same 
imagined viscosity. Qgilvie (119991) introduced a more complicated version of this approach 
by creating a quasi-linear theory of this model, one in which the angular mo men tum orien¬ 
tation diffusion coefficient became a function of several parameters. N els on & Papaloizou 
( 120001 ) constructed a formalism for relating the rate of radial transport of misaligned angu¬ 
lar momentum to the location of the orientation transition, but in the absence of a clear 
understanding about how the radial flows were regulated, t his f orma lism could not yield a 
definite statement about its location. Recently, iNixon & King (2012 ). also working within 
the isotropic viscosity formulation, pointed out that large misalignment angles might lead 
to an actual rupture of the disk, rather th an a s mooth transition. Shortly after, detailed 
hydrodynamics simulations (ISorathia et al.l 2013bf ) showed that the relation between angular 
momentum flux and disk warp has clear departures from the diffusion model. 

This situation can now be clarified. There is a well-established physical mechanism for 
internal stresses in accretion disks—correlated magneto hy drody namic (M HD) turbulence, 
driven by the magneto-rotational instability (iBalbus fe Hawlevi Il99ll 119981 ) . In addition, 
although treating MHD turbulence in disks requires very demanding numerical simulations, 
computer algorithms and computer power are now at a level such that even warped precess- 
ing disks can be treated directly. Disks tilted with respect to the rotation axis o f a centr al 


black hole have been simulated in fully relativistic three dimensional MHD (iFragile et al. 
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20071; McKinney et al.l 120131 Morales Teixeira et al.l 120141) . but the computational expense 
of full general relativity forced parameter choices (disk thickness, radial extent) that made 
accretion-related alignment processes difficult to see. The central difficulty is that if a nu¬ 
merical simulation is to follow the MHD turbulence, it must have a timestep very short 
compared to an orbital timescale, whereas the precession timescale where the orientation 
transition may occur is almost certainly many orbital periods long. Consequently, simula¬ 
tions able to probe the orientation transition at a realistic scale are prohibitively costly. 

However, it is possible to avoid this roadblock by using semi-Newtonian methods, in 
which the only relativistic effect retained is the lowest-order post-Newtonian term describing 
Lense-Thirring torque. Progress is possible in this framework because Newtonian gravity is 
both scale-free and obeys a linear field equation. As a result, in this treatment there is no 
maximum value for the spin parameter. Simulators are then free to choose the precession 
frequency associated with the Lense-Thirring torques, constrained only by the provisos that 
it must have the proper radial dependence, and it must be slower than the orbital frequency. 


By adopting this device, ISorathi a et al.l (j2013al ) were able to simulate a disk in which 
an alignment front propagated steadily outward. By studying this model in detail, they 
were able to achieve several things. First, they demonstrated that the internal MHD stresses 
triggered by disk warp are fundamentally not viscous in c haracter and h ave a magnitu de much 
smaller than predicted by the “isotropic a” model (M o rales Teixeira et al.l (2 0141 ) reached 
a similar conclusion). Second, although the magnetic forces associated with these stresses 
are generally weaker than the fluid pressure forces, the turbulence created by MHD effects 
significantly influences disk dynamics because it disrupts the propagation of bending waves. 
Third, they identified the specific mechanism by which alignment occurs: when the precession 
angle of the disk decreases outward (as it in general should for Lense-Thirring torques), the 
torque has a component that is in part opposed to the misaligned angular momentum at 
larger radius than where it is exerted. If the new angular momentum delivered to the disk is 
then transported outward by radial fluid motion, it progressively cancels the misalignment. 
Fourth, based on this understanding of the alignment mechanism, they developed a simple 
model for the speed of propagation of the alignment front that fit their data very well. 
However, because the radial surface density profile chosen for the simulated disk had a peak 
in the middle of the disk, there was too little misaligned angular momentum stored in its outer 
region s to be ab l e to slow, much less stop, the outward propagation of the alignment front. 
Thus, ISorath i a et ah] (12013a! ) were unable to study a time-steady orientation transition. 


In this pa per, we present a new simulation, very similar in most respects to that of 
Sorathia et al. (l2013al) . but with a surface density distribution altered so that the majority 
of the disk’s mass is located at large radius. This choice creates a situation in which the 
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alignment front can travel outward, but then decelerate and stop. This is just what occurs, 
and in this paper we report in detail on the properties of such a time-steady orientation 
transition. 


Simulation Details 


As in ISorathia et al.l (2013afh we use our Fortran-95 version of the 3D finite-difference 
MHD code Zeus flStone fe Norman!ll992al Jbll. The Zeus code solves the standard equations 
of Newtonian magnetohydrodynamics using direct finite differencing. The magnetic field is 
updated using the method of characteristics co nstra in ed transport (M OCCT) algorithm to 
maintain zero divergence to machine accuracy (IHawlev fe St oriel 119951 ). 


Our procedure is to begin with a conventional flat disk, evolve it long enough to allow 
the MRI-driven MHD turbulence to saturate, and only then turn on a torque-driving Lense- 
Thirring precession. To represent this torque, we add a force term of the form pv x h to 
represent the gravitomagnetic force per unit mass, where p is the mass density, v is the fluid 
velocity, and 

9.T f \(.T • 

(i) 


23 6(J ■ r)r 
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Here J represents the magnitude and direction of the spin vector of the central mass and r 
is spherical radius. We choose the direction of the spin vector to be in the x-z plane, tilted 
12° (0.21 radians) from the z-axis in the x direction. Its magnitude is such as to cause a 
test-particle at r = 10, our fiducial radius, to precess at a frequency 1/15 of the orbital 
frequency at that radius. 

The simulation grid uses spherical coordinates with a spatial domain in (r, 6, 0) that 
spans [4,40] in radius, [0.1, 0.9]7r in 6 and a full 2n in 0. Since we are using Newtonian gravity 
rather than (say) a pseudo-Newtonian gravity with a special lengthscale corresponding to 
th e gravitation al radiu s r g = GM/c 2 , the units of length are arbitrary. For consistency with 
Sorathia et ahl (l2013af) we use the same units, in particular measuring time in units of the 
orbital period at r = 10, and setting GM = 1. We use (352, 384,1024) grid cells in the radial, 
polar, and azimuthal directions respectively. The radial mesh is spaced logarithmically with 
a constant Ar/r_= 0.00654. We employ a polynomial spacing in the polar dimension (Eqn. 6 
of No ble et al. 02010), with £ = 0.65 and n = 13). This polynomial spacing focuses cells near 
the equatorial plane, placing 64% of the 6 zones within ±10° of the midplane. The </> zones 
are uniform in size. Outflow boundary conditions are employed on the radial inner and outer 
boundary, and along the 6 boundary that forms a “cut-out” around the polar axis. 


In our initial condition, the fluid orbits a point-mass in Newtonian gravity with a Ke- 
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plerian angular velocity, fl 2 = GM/r 3 . A similar initial disk was used in 
Its equation of state is isothermal, with sound-speed c 2 = 0.001, which corresponds to 0.1 of 
the orbital velocity at the fiducial radius of r = 10. The scale height of the disk, H = c s /f2, 
is then (r/10) 3 / 2 . If the scale height is instead defined in terms of a vertical density moment, 
the disk aspect ratio h/r ~ 0.07 at r ~ 6 and rises to ~ 0.12 at r ~ 20. The disk’s opening 
angle is thus about half the initial tilt of the disk relative to the black hole spin. We set 
the density p c — 1 at the equator at all radii, and the vertical distribution of density is 
P — p c exp(— z 2 /2H 2 ). The surface density £ increases with radius from the disk’s inner 
radius r; n = 6 out to r = 20 and declines outward from there to its outer radius r out = 28 
(see Fig. |T] for £(r) when the torque is turned on). At the inner and outer limits the disk 
is simply truncated. The initial condition is not in radial pressure equilibrium (especially at 
the disk boundaries). 

The initial magnetic field in the disk is a large set of nested dipole loops extending from 
the inner to the outer portions of the disk. The field is defined by a vector potential 

Ap = A 0 p 1/2 sm(2Trr/XH)(r/r in - 1)(1 -r/r out ). (2) 


Shi et al. (2012). 


The scalefactor A is set equal to 2 jc s . The vector potential is limited to positive values with 
a cutoff at 0.05p c , i.e., 

A ^ = rnax(A 0 - 0.05p c , 0). (3) 

The field amplitude factor A 0 is chosen so that the volume-integrated ratio of gas to mag¬ 
netic pressure, the plasma (3, is 1000. Turbulence is seeded by imposing random pressure 
perturbations at the 1% level on the initial condition. 


Throughout the simulation we monitor a number of quantities in order to gauge the nu¬ 


merical qual ity of the simul ation. These m e trics were developed and studied in lHawlev et al. 


(120111. 201311 . as well as in iSorathia et al.l (20121). One set of metrics measures the ratio of 
the characteristic MRI wavelength A mri to the grid resolution. Specifically, 


and 
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The characteristic MRI wavelength is defined as A mri = c 2' k \va\/Q, where the Alfven speed 
Va is obtained from the appropriate component of the magnetic field ( B g for Q e , B $ for Q^). 
The Q num bers are then t he nu mber of grid cells that sp an the characteristic wavelength 
of the MRI. Hawlev et ah ( 2011 1 and Hawley et. al. J2013 1 estimate that Q values > 15-20 
are indicative of adequate resolution. For the initial magnetic field given here, A mri ~ 0.1, 
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and the initial Qe values at the disk equator inside the fiducial radius range from ~ 4-8. 
There is no initial toroidal field, but orbital shear creates toroidal field from radial field on 
the orbital timescale. Thus, within a single orbit we can expect ~ B r , producing an 
azimuthal quality factor Q$ = 2-15 at r = 10. Although this would seem like marginal 
resolution, the Q values increase as B grows beyond the initial field amplitude (see below). 


A second set of metrics measures the average properties of fully developed MHD tur¬ 
bulence; these are, by definition, irrelevant in the initial condition of the simulation, so we 
do not evaluate them until the turbulence develops. These metri cs are calibrated by the 
values measured in highly-resolved local simulations. Hawlev et al . 120111. 2013!) develop two 
such diagnostics, a mag = M r ^jP mag ^ the ratio of the Maxwell stress M r $ to the magnetic 
pressure, and (Bp/(Bp, the ratio of the radial to toroidal magnetic energy. When suit¬ 
ably averaged over the computational domain in well- resolved s imul ation s, these quantities 
approach limiting values of 0.45 and 0.2, respectively flHawlev et al.l 201 31). 


3. Results 

3.1. Initial disk evolution 


The aim of our new simulation is to follow the process by which an accretion disk sub¬ 
jected to Lense-Thirring torques reaches a steady-state in which its inner portions are aligned 
with the central black hole while its outer portions remain oblique. Because the^ c haracter 


of this process depends upon the fluctuations created by MHD turbulence flSorathia et al. 


2013al) . we begin by evolving the initial disk until it is turbulent throughout most of its 


radial extent. Only after that state is achieved do we switch on the torque. 


Because the initial disk was not in radial pressure equilibrium, a period of adjustment 
occurs as the disk evolves. Although the unbalanced pressure gradients are small within 
the body of the disk due to the small ratio of sound speed to orbital speed, they are large 
at the disk boundaries. Initially, the inner edge of the disk pushes inward, but it then 
rebounds, driving a strong compression wave through the disk, while the outer radius of the 
disk expands outward, sending a rarefaction wave into the inner disk. 


By t ~ 12 fiducial orbits, which marks the end of the torque-free evolution, turbulence 
has been well established in most of the disk. Inside r = 20 at this time, the mass-weighted 
mean value of Qq varies between 15 and 22, Q$ is between 35 and 50, a mag has an average 
value of 0.3 7, and ( Bp /( Bp is 0.16 . The two Q values are well above the resolution thresh¬ 
olds set bv lHawlev et al.l ( 2011 . 2013|), while the a mag and (Bp/(Bp values are a little bit 
lower than the limiting values found in that paper. 
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Under the influence of the magnetic stresses associated with the turbulence, the disk 
spreads radially, both inward and outward. At the end of the torque-free evolution period, 
the disk has 65% of its original mass, and a surface density with the form sho wn in Figure [Q 
Unlike the radial surface density profile of the disk studied by lSorath i a et al.l (]2013aj), which 
reached a peak at r ~ 12 and declined sharply from there outward, the surface density in 
this simulation rises all the way to r ~ 20, reaching a peak roughly twice the surface density 
prevailing for 6 < r < 10. Thus, more than enough “misaligned” angular momentum 
Lj_ = L — L ■ s, with s the direction of the nominal black hole spin, is present on the grid to 
hold back the outward progress of an alignment front for a long time. 


3.2. Evolution with torque 

At 12.4 fiducial orbits, we turned on the torque and then ran the simulation for another 
12.4 fiducial orbits. Because the black hole spin implied by the torque is tilted 12° from the 
polar axis of the coordinates, the disk acquires an effective misalignment angle /3 = 12° at this 
point in the simulation; at subsequent times, we define the radius-dependent misalignment 
angle j3(r) = arcsin [L±(r)/L tot (r)], where L tot (r) is the magnitude of the shell-integrated 
angular momentum at radius r, and L±(r) is the magnitude of its component in the plane 
perpendicular to the black hole spin. 

If orbiting material moved as independent test-particles, the inner rings would precess 
much faster than the outer rings because the precession frequency declines outward oc r -3 . 
To measure the actual precession, we define the precession angle </> prec = tan ~ l (L' y /L' x ), 
where L' x are the shell-integrated x and y components of the angular momentum relative 
to a coordinate system found by rotating the grid coordinates around the y axis in order to 
make the black hole spin define the z axis. As Figure [2] shows, although the inner rings do 
precess faster, the radial gradient is noticeably shallower than r~ 3 : after 1 orbit, although a 
test-particle at r = 5 would have precessed by ~ 1.17T, the actual fluid ring at that radius 
has precessed by only ~ 0.67T, while the ring at r = 10 has precessed by an amount similar 
to the test-particle prediction; after 3 orbits, precession is ill-defined at r = 5 because that 
portion of the disk has already almost completely aligned, but the difference in precession 
angle between r = 10 and r = 15 is only ~ 0.17T, in contrast to the test-particle expectation 
Of ~ 0.37T. 

Several other features in this diagram are even more contrary to test-particle expecta¬ 
tions. Most importantly, the precession angle hardly changes from t > 15 onward for r < 10. 
Also beginning at f ~ 15, or after only ~ 2.5 orbits of torque, a zone at larger radius, 
7 1$ r ~ 13, reverses its sense of precession for several orbits. Once the sense of precession 
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Fig. 1.— Surface density as a function of radius at t = 12.4, the beginning of the torque. 
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at large radius resumes its normal sense (at t ~ 20), the precession phase gradient becomes 
very shallow throughout the disk. 

Meanwhile, the innermost disk aligns very swiftly, within a few fiducial orbits. The 
location of the alignment front moves rapidly outward, with the half-alignment point reaching 
r ~ 10 only ~ 3.5 fiducial orbits after the torques begin. However, there is a turn-around in 
the process and the front retreats over the next several orbits until it reaches its steady-state 
position at r ~ 7-8. From f ~ 20 until the end of the simulation at f ~ 25, neither the 
position nor the structure of the front changes appreciably. 


Even though the disk remains near its original alignment at r > 12 after t ~ 20, it 
continues to precess slowly at large radii, turning at a rate ~ 0.037T radians per fiducial orbit, 
roughly the test-particle precession rate at r ~ 16. Although there is little torque inside the 
alignment radius (r ~ 7), there is sufficient obliquity at larger radii that the disk continues 
to feel a torqurn This t orque is spread to larger radii by the same mechanisms seen in earlier 
studies f Sorathia et al. 2013bl ah When the radial precession gradient creates an inter-ring 
warp of nonlinear amplitude, it triggers a warp pulse that travels outward. Nonlinear in 
this context means that the parameter 0 = \M/d\nr\/(h/r) > 1, where £(r) is the unit 
vector in the direction of the total angular momentum associated with the shell at radius 
r. When that criterion is satisfied, the warp creates an order-unity radial pressure contrast. 
The resulting transonic fluid motions mix angular momentum between rings, smoothing the 
warp. This same smoothing process keeps the pulses narrow. Two groups of such pulses 
can be clearly seen in Figure SJ During the first ~ 3 orbits, five pulses of strong (-0 ~ 4-6) 
warp propagate outward. Although each one travels at very nearly constant speed, the pulse 
speed declines from first to last, falling from ~ 8 to ~ 2, i.e., from ~ 1.3c s to ~ 0.3c s , from 
first pulse to last. Their direct impact on local precession angle can be seen immediately by 
comparing Figures [2] and HI A second group, somewhat weaker (0 ~ 2-3), are visible in the 
zone 8 < r < 15 and 20 < t < 25. Although less dramatic, their effect on precession phase 
can likewise be seen by examination of the region in radius-time space occupied by these 
pulses in Figure [2J All of this group of pulses have rather similar speeds, ~ 3 length units 
per fiducial orbital period, ~ 0.5c s . As hinted by this speed, these pulses are related to linear 
bend in g waves, w hose characteristic propagation speed is half the isothermal sound speed 
flLubow et al. 2002). However, a number of factors combine to make their behavior somewhat 
different. In addition to the rapid local damping caused by their nonlinearity, Figure [4] also 
shows that their amplitude is quite irregular, a symptom of their propagation through a 


backgr ound medium made turbulent by the magneto-rotational instability flSorathia et ah 


2013al h 


Throughout the evolution, the disk continues to accrete due to internal stresses, but in 
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Fig. 2.— Color contours (see color bar) of the precession angle (p prec /7T as a function of radius 
and time. Although the simulation extended to r = 40, we show only the range 4 < r < 25 
because little happens at larger radius. The precession angle in the small region at r < 7 
and t ~ 14 is ill-defined because the disk is almost exactly aligned (see Fig. [3]). 
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Fig. 3.— Color contours (see color bar) of the misalignment angle /3/7T as a function of radius 
and time. Although the simulation extended to r = 40, we show only the range 4 < r < 25 
because little happens at larger radius. 
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Fig. 4.— Color contours (see color bar) of the normalized disk warp i/j. 
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the vicinity of the alignment front, the radial motions induced by the disk warp create a 
large enough Reynolds stress to enhance the local inflow rate. To measure this effect, we 
compute the mean radial, i.e., accretion, velocity from the mass flux as follows: 

, r \ _ I /wV 2 sin ddOdcf) 
f pr 2 sin ddddcj) 

This velocity is then averaged over time for the last orbit of the simulation. The result is 
shown in Figure 0 as the ratio of the accretion velocity to local orbital velocity. This ratio 
is ~ —0.004 well outside the transition front (13 < r < 18), but its magnitude increases 
inward from that radius, becoming ~ 50% greater by r ~ 7, the location of the front. 
Another consequence of this effect is a shift in mass from the surface density peak, which is 
at r ~ 15-20 when the torque begins, to the transition front region (compare Figs. [0] and [H). 
Note that the diminution in radial speed at very large radius is an artifact common to any 
accretion simulation with a finite disk: matter on the outside must move still farther away 
as it absorbs angular momentum transported to larger radius by the accretion stresses. 

Over the course of the evolution with torque, the quality parameters Qq^ maintain 
values consistent with adequate resolution. Within the region of the disk running from 
r = 6-20, Qe ranges from 9-14, while Q^ varies between 20 and 30. 


4. Analysis 
4.1. Overshoot 


Because the torque is turned on abruptly, the disk is subject to some transient effects 
as it adjusts; one of these is an overshoot in the location of the alignment radius. This 
“overshoot”—the quick progress outward of the alignment front followed by a retreat to 
its long-term position—can be understood as a combination of several previously-identified 
processes. 


The overshoot begins as the five strong warp pulses described in the previous section 
propagate outward. When the inner regions are advanced in precession phase relative to the 
outer regions, the gravitomagnetic torque delivered to the inner regions is dir ected so that it 


can pa rtially cancel the misaligned angular momentum of the outer regions flSorathia et al. 


2013al ). Accomplishing this cancellation requires transport outward, and the pulses do this. 


Alignment by this mechanism can be particularly rapid when there is little warp at still 
larger radii because radial mixing flows bringing inward angular momentum aligned with 
the original orientation of the disk are weak. 
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Fig. 5.— Ratio of radial velocity to orbital speed as a function of radius. Negative radial 
velocity indicates inward motion. 
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Fig. 6.— Surface density as a function of radius at t — 28, shortly before the end of the 
simulation. 
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However, a consequence of quick alignment is the creation of significant warp. After 
roughly a local orbital time (and at the relevant radii near r ~ 10 , the local orbital time is 
the sa me as our fiducial orbital period), the newly-created warp accelerates radial mixing 
flows ( Sorathia et ah 2013bl b The delay time is roughly a local orbital period because this is 
the local dynamical time. These radial flows carry inward angular momentum that is both 
more misaligned and further behind in precession phase. As a result, the local alignment is 
reduced and the precession temporarily reverses. Thus, both the alignment front overshoot 
and the precession reversal can be understood as consequences of transients associated with 
the beginning of torque in a flat, but misaligned, disk. Although it is somewhat artificial 
to turn on the Lense-Thirring torque abruptly, this procedure does allow us to highlight the 
dynamic alignment process as it subsequently occurs. The disk ultimately achieves a steady 
state alignment configuration, despite the impulsive initial nature of the torque. 


4.2. Location of the steady-state orientation transition 


The ultimate location o 


the a lignment front is also very consistent with the formalism 


proposed by ISorathia et al. Q2013aJ) to estimate both the propagation speed of an alignment 
front and the location where it stalls, reaching a steady-state. They proposed that the 
alignment front should move outward at a rate 

drf 2(cos 7 )a*(GM) 2 


, rf 


dt 


sin f3(rf)c 3 r 3 J 2 Yj(rf) Jo 


dr' sin /3ir') E(r')/ r ' 3 ^ 2 , 


(7) 


where 7 is the mean angle between the angular momentum being mixed outward and the 
local misaligned angular momentum, and f3 is the local misalignment angle. Testing this 
model against their simulation data, they found excellent agreement. The expression for the 
front speed can also be rewritten in an entirely equivalent form, in which the integration 
variable is transformed to the dimensionless x = r'/rf: 


drt . . , f 1 _ 3/2 sm/3(x) E(x) 

— = (cos 7 )r/O prec r/ / dxx ' 7 . ■ 

dt Jq smp(r) E(r) 

In other words, the characteristic speed of the alignment front is rs./ r f Q 


( 8 ) 


prec* 


Sorath i a et al. (l2Q13af ) suggested that the place where the front should stall is defined 


by the place at which the speed of the alignment front is matched by the speed with which 
misaligned angular momentum from the outer disk can be mixed inward. From their simula¬ 
tion data, it was impossible to estimate the speed of inward mixing. However, they proposed 
a simple parameterization for describing it, in which the speed of misaligned angular mo¬ 
mentum flow was taken to be ry/t in . In this language, the steady-state orientation transition 
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would be found at 77 = i?r, defined by 

^prec(-Rr)bn = ((cOS7)X)^ 1 , (9) 

where X is the dimensionless integral in equation [BJ 

Our simulation data give the location of Rt, so equat ion [9] can be inverted to find the 
characteristic inflow time for misaligned angular momentum: 

tin = ((cos 7 )XO prec (Rt ))- 1 . (10) 


Taking Rt — 7, we find that during the period 18 < t < 25, when the alignment front 
appears to be almost constant in position, t- m ~ 11.5 fiducial orbits, the exact value varying 
slightly from time to time. As a standard of comparison, we note that the orbital period at 
r = 7 is ~ 0.6 fiducial orbits. 


One possible in terpretati on of this mixing time is in terms of a diffusion model (jPapaloizou fe Pringle 


1983l : Pringle 119921). although iSorat hia et al. fl2013bl ) showed that this is not a perfect match 
with actual time-dependent response to disk warping. However, even if it does not de¬ 
scribe time-dependent effects, a diffusion model might be an appropriate device with which 
to model the radial angular momentum mixing supporting a time-steady alignment front. 
If one does describe the mixing as a diffusion process, the associated diffusion coeffici ent 
should be ~ c 2 /D because the characteristic speed of radial flows v r ~ c s (ISorathia et al. 


2013bl) . and orbital mechanics limits the distance such a flow can travel radially to 
v r /n h. Although dish 


bv PaDaloizou & Prina 

e 

(1983(1. the associated radial velocities are not regulated bv an 

“isotropic a viscosity” ( 

Sorathia et al. 2013a 

Morales Teixeira et al. 2014b but bv hvdrodv- 

namic expansion, gravity, and weak shocks ( 

Sorathia et al. 

2013b 

). The characteristic speed 


of mixing is then ~ $(h/r) 2 (r/Ar)rQ(r), where A r/r is the characteristic lengthscale for 
orientation change within the front, and $ is a number of order unity that can be determined 
from simulations like ours. Defining Ar = \A£\/\d£/dr\ for A£ the difference in direction 
between the disk’s initial orientation and its aligned direction, we find Ar ~ 6. Equating 
this speed to the unimpeded drj/dt leads to an estimate for the stationary front location 


Rt — R* 


/ \ <-7-^prec(-^*) At/T 

(cosq/X- 


fi(i2*) $(h/r) 


2/3 


( 11 ) 


In order to find a scale for Rt, we have created uni ts relat i ve to a f iducial radial scale R *. 
Note that in terms of the formalism created by iNelson & Papaloizou (1200011 . this expression 
corresponds approximately to the geometric mean of Tsfi and 7 sf2- 
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If this model correctly describes the dynamics of alignment front equilibrium, the order 
unity quantity $ is 


$ ~ (cos7)X[O prec ( J R T )/0( J R T )] (Ar/r)(h/r) k (12) 

Our data indicate that <f> ~ 0.6-0.8, depending on the time at which it is measured. 

It is possible that $ might be independent of other parameters; however, that remains 
to be demonstrated. It could be, for example, that $ is dependent on c s /v or b ~ H/r. It 
should also be noted that because the width of the orientation transition was smaller than 
the disk simulated, the transition thickness Ar and the internal warp rate i/j were free to 
adjust to whatever the dynamics required. There may be cases in which the intrinsic disk 
obliquity is large enough, and the radius at which the disk is fed is small enough, that 
additional constraints are placed on Ar and therefore 0; these might then also alter $. In 
extreme cases, shock heating might be strong enough to alter c s , leading to nonlinearity in 
the effective diffusion rate. 


4.3. What is the role of bending waves? 


As discussed in Sec. 13.21 we see clear evidence for nonlinear bending wave pulses moving 
through the aligning disk. However, these pulses propagate in a rather irregular way, their 
amplit udes fluctua tin g stron gly as they move outward through the disk. In our previous 
paper (jSorathia et al. 2013al ). we showed that in a purely hydrodynamic situation these 
pulses propagate much more regularly. We attributed the contrast to the fact that warped 
hydrodynamic disks are laminar, whereas warped MHD disks are turbulent; the turbulence 
disrupts, but does entirely destroy, these waves. One consequence (as seen both in the 
previous simulation and in this one) is that an MHD disk is not able to maintain a state of 


solid-body rotation, which it does achieve in a purely hydrodynamic state (ISorathia et al. 
2013at) . 


Since the work of iPapaloizou Sz Pringle (119831) . it has been customary to categorize 
the dynamics of warped disks in terms of a division between a “bending wave” regime, in 
which the ratio of vertically-integrated accretion stress to vertically-integrated pressure a is 
< H/r , and a “diffusive” regime, in which a > H/r. It might then seem a natural question 
to ask how the phenomena we observe may be classified in these terms. Before doing so, 
however, it is useful to recall the origin of this distinction. It was derived by assuming that 
linear bending waves damp at a rate af2, while the speeds of the radial motions induced by 
the warp-created radial pressure gradients are controlled by an “isotropic viscosity” acting 
to reduce vertical shear. “Isotropic” meant that if the vertical shear were the same as the 
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orbital shear, the r-z component of the stress would be the same as the accretion stress, the 
r-<f) component. 

However, neither of these assumptions fits well with the conditions in warped MHD 
disks. In the vicinity of an alignment front, the warp is almost always nonlinear in the sense 
that ip > 1 because even the transonic radial motions induced by nonlinear warps cannot 
spread the front enough to make it linear. This condition is generic because the location of the 
front is determined by the condition that the mixing timescale is comparable to the precession 
timescale. Indeed, the only regions in our simulation in which the warp is linear are at 
sufficiently large radius that they are essentially untouched by the torques. Co nseq uently, the 


relevan t damping rate is the one applicable to nonlinear waves. As shown in ISorathia et ah 


(J2013b|), damping of nonlinear warps takes place very rapidly, and the mechanism is exactly 
what is commonly associated with the “diffusive” regime: radial motions driven by the 
radial pressure gradients associated with the warp mix differently-aligned angular momentum 
across a radial orientation gradient. Thus, nonlinear bending waves combine aspects of these 
two, nominally contrasting, regimes. 


In addition, as shown in ISorathia et al .1 (I2013al) and Morales Teixeira et al.l (1201411 . there 
is no “isotropic viscosity”. The r-z component of the Maxwell stress does not behave like a 
viscosity: it is as likely to strengthen shear as to weaken it. Moreover, its magnitude, when 
measured in pressure units, is typically three orders of magnitude smaller than the stress 
associated with accretion. Consequently, the speeds of the radial motions are limited not 
by a “viscosity”, but, when the warp is nonlinear, by the dynamics of free expansion in the 
presence of gravity and by weak shocks that occur when the radial motions encounter other 
parts of the warped disk. For this reason, the radial speeds in regions of nonlin ear warp 


are in general transonic, with Mach numbers that increase with greater ip (ISorathia et al. 
2013bj). Thus, it is questionable whether the “bending wave” vs. “diffusive” distinction can 


be applied here; the only distinction having the potential to be qualitative is between linear 
and nonlinear bending waves. 


4.4. Do warped disks break? 


When disks subjected to Lense-Thirring torques are simulated in pure hydrodynamics 
with an assumed isotropic viscosity and using the SPH algorithm, under certain circum¬ 
stances breaks have appeared between the outer misaligned secti on and the inner aligned 
portion (INixon fe Kindl2012l: iNixon et al .1120121: iNealon et al.ll2015l) . One proposed criterion 
for this to happen is that the external torque 7~ exceed th e inter-ring torque associated with 
accretion stresses G (INixon et al .1120121: INealon et al.ll2015l) . Alternatively, it may be that the 
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cr iter ion is for theprecession timescale to be shorter than the radial sound-crossing timescale 
f Nealon et al. 2015 ). 


We use several diagnostics to test our late-time data for disk breaks in this simulation. 
The first is continuity in the azimuthally-averaged surface density. As can readily be seen 
in Figure |6l £(r) varies quite smoothly with radius. Indeed, the surface density near the 
alignment transition hardly changes over the ~ 15 fiducial orbits in which the torques op¬ 
erate, while the surface density at somewhat larger radii increases. Another is the density 
distribution in the poloidal plane measured at the azimuth corresponding to the outer disk 
precession angle. This is shown in Figure [7] The plane of the disk does vary, of course. At 
small radii, where it is aligned normal to the black hole spin, it is tilted by 12° with respect 
to the grid; at large radii, precession around the black hole spin direction moves it out of the 
equatorial plane. Throughout the radial extent of the disk, the density distribution is lumpy, 
but the maximum contrast in density within the midplane is only about a factor of 2. Thus, 
this view also provides no evidence of breaks. A third diagnostic is the density distribution 
on a spherical shell at the radius where the warp is greatest (Fig. [8j). Because the mean an¬ 
gular momentum on this shell is tilted by a bit more than half the initial inclination relative 
to the coordinate equatorial plane, and it has precessed, the disk plane traces a sinusoidal 
pattern in polar angle relative to the 6 x 0 plane in this coordinate system. It shows density 
fluctuations at the factor of 3 level, but there is also no indication of a dramatic break here 
either, even though this is the radius where ^ is greatest (see Fig. [4]). 

It is therefore instructive to test whether the proposed break criteria apply to a calcu¬ 
lation in which all stresses are physical and numerical transport can be readily controlled. 
The first criterion can be concisely phrased if the actual accretion stress is measured in 
pressure units, i.e., in terms of a, the time-averaged ratio of vertically-integrated stress to 
vertically-integrated pressure. It is then 

T (O prec /Q) sin 6» 

G ~ a(h/r) 2 ' 1 j 

At the location of the orientation transition (r ~ 7), the frequency ratio is ~ 0.1, while 
our measured stress ratio a declines with radius from ~ 0.1 at r ~ 5 to ~ 0.05 at r ~ 10. 
Combined with our disk aspect ratio h/r ~ 0.1 and the outer disk inclination sind ~ 0.2, 
we find T/G ~ 27. Thus, even though the precessional torque is considerably greater than 
the inter-ring accretion torque, the disk does not break. 

One possible explanation for this discrepancy is that the argument for the T/G criterion 
assumes that ordinary accretion stresses are responsible for restoring matter to regions in 
which the precessional stress might otherwise create gaps. However, warped disks induce 
strong Reynolds stresses, and these can easily be much stronger than the Maxwell stresses 
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Fig. 7.— Linear color contours (see color bar) of the disk density in a poloidal slice at 
(j) = 0.37T, the approximate precession angle of the disk at t — 24. 






theta 
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Fig. 8.— Linear color contours (see color bar) of the disk density at t = 24 on the radial 
shell at r = 7.5, the radius of sharpest warp. 
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due to correlated MHD turbulence that ordinarily drive accretion. In fact, at a steady-state 
alignment transition front, these Reynolds stresses produce a torque that is—by definition— 
equal to the aligning torque at the front, and the aligning torque is comparable to the full 
Lense-Thirring precessional torque. Thus, if the correct criterion for producing a break is 
for the Lense-Thirring torques to exceed the smoothing torque due to internal stresses, there 
could never be a strong break at a disk’s most vulnerable point, the alignment transition. 

The second criterion amounts to computing the product 

^prec/^ ( ' 14 '| 

a L prec Across 7 / ? V 

n/r 


where t CTOSS is the time for sound waves to traverse a distance comparable to the radius where 
this product is measured. If the theory for the location of the transition radius outlined in 
the previous section is correct, this criterion may be rewritten as 


n 


prectcross 


$ 

(cos7)X(Ar/r) 


(15) 


Thus, this ratio is always ~ Z . For our simulation’s parameters, it is ~ 1 at the alignment 
transition. Just as for the previous criterion, our theory predicts that this criterion, too, 
should always be marginal, and therefore strong breaks would not be expected. 


Nealo n et al. (120151 ) suggest that this criterion applies in what they deem the “bend¬ 


ing wave” regime, a regime in which the pressure-coupling of bending waves provides the 
strongest coupling between adjacent rings. Because these waves generically travel at 0.5c s , 
the sound wave crossing time is a reasonable estimate of their crossing time also. As we 
have already argued, linear bending waves do not appear to be particularly important con¬ 
tributors to tilted disk dynamics. However, nonlinear warps result in transonic flows, whose 
effects also travel at ~ c s . In this sense, it is not surprising that we would predict that 
this criterion, like our reworked version of the first, should also always be marginal at the 
orientation transition radius. 


4.5. Is apsidal precession relevant? 



perhaps affecting overall disk alignment. Nonetheless, in our lowest-order PN accounting 
of general relativistic effects, we have included the gravitomagnetic term, but not the term 
associated with apsidal precession, even though the two effects are formally the same order in 
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the expansion. We did so for a number of reasons. One reason is that there is no way to make 
the apsidal precession consistent with the magnitude of our amplified gravitomagnetic term. 
However, we also believe it to have at most minor physical influence. Several arguments 
point in this direction. First, the amplitude of the velocity perturbations created by this 
effect is very small. They are proportional to the eccentricity e of the orbits, and e C-s / ^orb- 
As a result, the Mach number of these motions is only of order the ratio of the gravitational 
radius r g to the radius, which is generally quite small for r ~ R ?. They are then far slower 
than the transonic radial motions themselves; in fact, they are very small even compared 
to the fluid motions associated with the turbulence, which are ~ 0.1c s . Second, the radial 
flows are altered by hydrodynamics on the orbital period, as streams encounter other por¬ 
tions of the disk, sometimes through shocks. By contrast, during one orbital period, the 
magnitude of the apsidal precession angle is only 6nr g /r <C 1 for e -C 1. Thus, during the 
time over which a specific elliptical orbit exists, it is rotated by only a very small amount, 
while that fluid’s motions are changed at the order-unity level by other mechanisms. Third, 
the apsidal precession can also be viewed as due to a small offset between the frequency 
of radial epicyclic motions and the orbital frequency. However, in an MHD orbiting fluid, 
the magnetorotational instability makes radial epicyclic motions unstable, and the fastest 
growing mode has a growth rate ~ Q. It is hard to see how the apsidal precession would be 
significant relative to all these other effects. 


5. Conclusions 


Our simulation has shown that an accretion disk whose intrinsic orbital plane is inclined 
obliquely to the equatorial plane of the central spinning black hole can achieve a steady-state 
in which the inner disk aligns with the black hole’s equatorial plane while the outer disk 
retains its original orientation. 


Using the alignment front propagation formalism developed in a previous paper flSorat hia et al 


2013aJ), we have, for the first time, measured directly the timescale for misaligned angular 


momentum mixing due to the radial flows induced by disk warping. In this specific case, 
that time scale is ~ 12 local orbital periods. In the context of a time-steady alignment front, 
it is possibl e to estimate this mixing tim esca le in terms of a diffusion model of the sort first 
explored by Pa palo i zou fe Pri ngl e (119831) and Pringle! (Il992f ). Because nonlinear warps gener- 
ically produce transonic radial motions, one would expect the diffusion coefficient in such a 
model to be <hc 2 /f2, where $ is a dimensionless number of order unity. Interpreting the mixing 
time we measure in terms of such a model, one would then expect fr 1 = < h(/i/r) 2 (Ar/r)~ 1 fl 
for fractional front width A rjr. Calibrating this relation with our simulation data, we fold 
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that ~ 0.6-0.8 for our measured fractional front width A r/r ~ 6/7. Finding a coefficient 
of order unity is consistent with the thought that a diffusion model of this character may be 
appropriate to estimating the equilibrium properties of disk warps, even though this model 
is inadequate to describe their time-dependent properties. Stronger evidence for this propo¬ 
sition must await simulations with different values of the sound speed to orbital speed ratio 
and perhaps also different disk inclination angles; it remains possible that <f> could depend 
on these, or other, parameters. 


Nonetheless, if this estimate is good to factors of several, it permits estimation of steady- 
state orientation transition radii, given the black hole mass and spin and the accretion disk 
surf ace dens ity a nd as pect ratio profiles. In fact, this approach has already been taken by 


Mil l er fc Krolik (1201311 . who assumed ~ 1; the results of our new simulation support their 


assumption. 


We have also probed several quantitative aspects of the alignment front’s approach to 
equilibrium. For example, we observed an overshoot in the front’s location before settling 
into a steady-state. A pure diffusion picture of warp motion would not yield this sort of 
behavior, but its qualitative character (in particular, its re laxation t ime) is cons istent with 
previous detailed studies of warped disk hydrodynamics (ISorathia et ah 2 013b ). We have 
also begun to constrain suggestion s that disks underg oi ng Lense-Thir r ing to rques can suffer 
sharp breaks (INixon fe Kind 120121; iNixon et al.l 120121) . iNealon et al.l (1201511 offered criteria 
for when these breaks should develop; despite satisfying one of their criteria strongly and 
another one marginally, we fold no evidence for this happening. If their primary criterion 
were reframed to include the effects of the radial motions induced by a strong warp, we 
would find both are marginal, and for the same reason: at a steady-state alignment front, 
the torque due to the Reynolds stress of these motions is automatically matched to the 
torque applied by gravitomagnetism, thereby creating a relatively rapid change in disk plane 
orientation across the front, but not a sharp break. 


A number of questions remain open for future work. Among the most interesting are 
the degree to which <h may vary with parameters such as c s /v or b and the possible effect of 
boundary conditions that prevent the transition front from extending across as wide a radial 
range as internal dynamics would demand. If the location of disk-feeding is at a relatively 
small radius and the total change in angle across the alignment front is relatively large, 
the warp rate ^ may be forced to become larger than it would be otherwise. In such a 
situation, there might be both a significant alteration in <f> and enough shock heating to 
create nonlinearity through thermodynamic effects. 
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